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Abstract. We explore the initial conditions for cosmological N-body simulations 
suitable for calculating the skewness and kurtosis of the density field. In general, the 
initial conditions based on the perturbation theory (PT) provide incorrect second-order 
and higher-order growth. These errors implied by the use of the perturbation theory to 
set up the initial conditions in N-body simulations are called transients. Unless these 
transients are completely suppressed compared with the dominant growing mode, we 
can not reproduce the correct evolution of cumulants with orders higher than two, even 
though there is no problem with the numerical scheme. We investigate the impact of 
transients on the observable statistical quantities by performing A^-body simulations 
with initial conditions based on Lagrangian perturbation theory (LPT). We show that 
the effects of transients on the kurtosis from the initial conditions, based on second- 
order Lagrangian perturbation theory (2LPT) have almost disappeared by z ~ 5, as 
long as the initial conditions are set at z > 30. This means that for practical purposes, 
the initial conditions based on 2LPT are accurate enough for numerical calculations of 
skewness and kurtosis. 
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1. Introduction 

The evolution of tlie large-scale structure in the Universe is one of the most important 
topic in astrophysics. The standard scenario for the structure formation is that the 
primordial density fluctuation grows through its gravitational instability [U El El HI |5] . 
Even though the perturbative descriptions are possible when the density fluctuation is 
small enough, in order to follow the distribution far into the nonlinear region, we must, 
inevitably, rely on the numerical calculations, named, iV-body simulations [6]. 

For A^-body simulations, there is a problem about how to set up the initial 
conditions. Even though the naive expectation is that it is better to start simulations as 
early as possible, like recombination era, it is well known that this is extremely difficult 
numerically. When we start A^-body simulations at an early era, the initial condition 
can be overwhelmed by sources of noise, such as the numerical roundoff error or the 
shot noise of the particles. Therefore, in the standard scheme [7], we must use the 
perturbative approach. Until the fluctuations come into the quasi-nonlinear regime and 
stable numerical calculation become possible. 

Even though there has been much progress about A^-body code with how to solve 
the non-linear structure in the universe, the method to set up the initial conditions 
into these codes has been almost the same from the first works of this field [El |9] . In 
many cases, the Zel'dovich approximation (ZA) [10], i.e., the first-order approximation 
of Lagrangian perturbation theory (LPT) have been applied for the initial conditions 
of A^-body simulations for a long time. (For reviews of LPT, see for example [HI [12]). 
Even though the ZA describe the growth of the density fiuctuation much better than the 
Eulerian linear theory, it does not reproduce the higher order statistics, like skewness 
and kurtosis with very poor accuracy [T3l [TH [15], because the acceleration is always 
parallel to the peculiar velocity. In other words, the acceleration does not refiect the 
exact clustering in the ZA. 

Therefore, A^-body simulations with the ZA initial conditions fail to pick up the 
correct second- and higher-order growing modes hence failing to reproduce the correct 
statistical properties of the density fiuctuation until very late times [HI [T7] . 

For this problem, recently, Crocce, Pueblas, and Scoccimarro [18] proposed the 
improvement by adopting different initial conditions. Basically, their initial conditions 
are based on the approximations valid up to second-order Lagrangian perturbation 
theory (2LPT) which reproduce exact value of the skewness in the weakly nonlinear 
region [191 120] • With these initial conditions, they calculate the statistical quantities 
and show the effects of transients related with 2LPT initial conditions decrease much 
faster than the ones related with ZA initial conditions, that is, the transients with 2LPT 
initial conditions are less harmful than ones with ZA initial conditions. 

However, there still exist transients related with 2LPT initial conditions which 
prevent to reproduce the exact value of higher order statistical quantities like the 
kurtosis, and there is no guarantee that 2LPT initial conditions are accurate enough for 
these quantities. 
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Therefore, as a natural extension of [I^, we examine the impact of transients from 
initial conditions based on 2LPT in A^-body simulation. In this paper, in addition to the 
ZA and 2LPT, we also calculate the non-Gaussianity with the initial conditions based 
on third-order Lagrangian perturbation theory (3LPT) which reproduce exact value of 
the kurtosis in the weakly nonlinear region [211 ES], [231 EH EH] . 

This paper is organized as follows. In Sec. [2l we present Lagrangian perturbative 
solutions valid up to the third-order and briefly explain the origin of transients. Then, 
after introducing the statistical quantities of interest in this paper in Sec. [31 we show 
the methods and the results of numerical simulations in Sec. [H For these results, we 
provide alternative interpretation based on simpler model in Sec. [5l Sec. [6] is devoted 
to conclusions. 



2. Lagrangian perturbations 

In this section, we briefly summarize Lagrangian perturbation theory (LPT) and obtain 
the approximations valid up to the third-order. We also point out explicitly why 
transients appear in the perturbative approach. 

Before considering the perturbation, we present the background cosmic expansion 
which determines the motion of cosmic fluid in Newtonian cosmology. We consider 
ACDM model in which the Friedmann equations are given as 

Id^a AtvG a , , 

with a background density ph of pressureless fluid and a cosmological constant A. 

Even though we take account of A for the numerical calculations, we sometimes 
consider the case with A = 0, so-called Einstein-de Sitter (E-dS) model, for a concrete 
solution in this section and the next section. This can be justified since the effect of A 
is negligible at the time we set initial conditions for cosmological A^- body simulations. 
For simplicity, in this paper, we do not consider back-reaction from the motion of the 
matter to cosmic expansion. 

Next, we consider the perturbations. In this paper, we consider the Lagrangian 
perturbation in which solutions for cosmic fluid are already derived by several people 
[TOl [T9l [201 l22l [231 El]- For this purpose, it is necessary to define the comoving 
Lagrangian coordinates q in terms of the comoving Eulerian coordinates x as: 

q = x + S{x,t), (3) 

where S is the displacement vector which denotes the deviation from the homogeneous 
distribution. 

It is worth noting that it is not the density contrast S but the displacement vector 
S that is regarded perturbative quantity in LPT. 
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In the Lagrangian coordinates, from the continuous equation, we can express the 
density contrast exactly as 

5 = J-i-l, (4) 

where J is the determinant of the Jacobian of the mapping between q and x: dx/dq. 
From Eq. (jlj), we can see that the break down of the perturbation with respect to 5 
does not necessarily mean that to which is the strong reason for considering the 
Lagrangian picture. In other words, although the density diverges or the fluid forms 
Zel'dovich pancake, the perturbation does not diverge. 

From the physical property, S can be decomposed to the longitudinal and the 
transverse modes: 

5 = 5^ + 5^ , (5) 
V, x5^ = 0, (6) 
V,-5^ = 0, (7) 

where Vg means the Lagrangian spatial derivative. 

In this paper, since it is well known that the transverse mode is negligible for 
pressureless fluid, we consider only longitudinal mode. From Kelvin circulation theorem, 
this is quite natural if it is generated only by the action of gravity. For the discussions 
with the transverse mode in LPT, see [26] . 

Hereafter we mainly follow the description by Catelan [21]. First, we consider the 
perturbations in Einstein-de Sitter Universe. In the Lagrangian description, changing 
the temporal variable as r = the basic equation for the density contrast for the 

pressureless fluid named Lagrangian Poisson equation is given by 

[(1 + V, ■ S)6^p - S^p + S%] Sp^ = a(r) [J{q, r) - 1] , 

(8) 

where Sap = dSa/dqp is the deformation tensor, S'^p is the cofactor matrix of Sap- Dots 
denote the differentiation with respect to r and a(r) is a function of r which includes 
the information of cosmic expansion law. Especially for E-dS universe, oi{t) = 6r~^. 

We now solve the dynamical equations for the displacements S according to the 
following Lagrangian perturbative prescription: 

Siq, r) = 5(1) (q, r) + S^^) {q, r) + S^'^ (g, r) + • • • . (9) 

Here S''-^) corresponds to the first-order approximation, S^'-^^ to the second-order 
approximation, and so on. Since we consider only the longitudinal modes, the 
perturbative solutions are described with gradient of scalar functions, 

5(") = V,^(") . (10) 

For the pressureless fluid, it can be shown that the perturbative solutions are 
separable into the temporal and the spatial parts, 

5«(g,r) = Z}(r)5«(q), (11) 

S^'\q,r) = E{r)s^'\q), (12) 

S^'\q,r) = F{r)s^'\q). (13) 
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The dynamics of the evolution constrains in general both the temporal dependence as 
described by the functions D, E, F, . . . and the spatial displacements s^"\ 

First we derive the Lagrangian linear perturbative solution (ZA). Making use of the 
fact that up to the third-order, the following expression for the Jacobian determinant J 
holds [2lj: 

1 



J = 1 + V, • 5(^) + ■ 5(2) + 



(V,-5«)^-5iXi 



+ V,-5(^)+[(V,.5«)(V,-5(^))-5lXi 
+ det(^«), 



(14) 



we can easily find the first-order approximation truncating Eq. ([8]), which yields, 

D{t) - a{T)D{r) = 0, (15) 

and no constraint to s^^\q). 

For the E-dS model, by substituting Q;(r) = 6r^^ into Eq. (IT^ . we can obtain the 
concrete form of the growing mode and the decaying mode solution as 



D^t) 



to) 

t 

To 



(16) 



(17) 



In general, we consider only the growing mode because the decaying mode is suppressed 
by a^^/^ compared to the growing mode and soon becomes negligible. 

It is worth noting that the analytic solution for the linear perturbative solution is 
obtained even in the case with A. The perturbative solutions are described as 



D_(t) = h, 



5 2 
6' 3 



where -B1//12 is incomplete Beta function: 



p, 



1 -pf-^dp 



(18) 
(19) 
(20) 

(21) 



Hereafter we consider growing mode (-D+) only and re-define D as D+. Even though the 
impact is not so significant, various modes appear in higher-order perturbation when 
we consider D_. For detail, see [27] . 

Next, we construct the solution valid up to second-order Lagrangian perturbation 
theory (2LPT). Retaining only the quadratic terms in Eq. ([8]) and using the first-order 
results, the system of equations become as follows: 

E{t) - a{T)E{T) = -a(r)D(r)2 
1 
2 



(V-s 



(1)^2 



(22) 
(23) 
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where s^^p = dsa/dqp. From Eqs. (122]) and we see that unhke the hnear case, the 
second-order approximation constraints both the temporal and the spatial dependence 
of the solution. 

In the E-dS model, by substituting a(r) = 6r~^ into Eq. (122!) . the growing mode 
solution can be obtained as 

E^{t) = -l('-X'' = -lD{tf. (24) 



7 v^o; 7 

In addition to this, the solutions satisfying Eq. fl22l) without the source term like 

E_{t) = ci^/^ + C2t-\ (25) 

where Ci and C2 are constants, also satisfy Eq. (1221) . 

Of course, E_{t) can be regarded as the "decaying modes" and after some time 
become negligible regardless of the concrete value of ci and C2 compared to the growing 
mode E^{t). They are, however, suppressed by only a~^, and if the initial conditions 
do not use the exact value for ci and C2, the error survives until late time. These 
nonphysical decaying modes are well known transients. 

Actually, if the initial conditions for A^-body simulations are set by the ZA, these 
are appropriate at only linear-order and are inappropriate at second- and higher-order. 
Therefore, we cannot obtain the correct ci and C2 as long as the initial conditions are 
set by the ZA. This is the point we discuss in this paper. (For the complete discussions 
about the origin of transients, see e.g. Sec. 2.1 in [18].) 

On the other hand, from Eq. (|23l) the spatial part can be described as 

s^''> = \ [s^'^ (V, ■ - ■ V,) 

+ , (26) 

where RP''' is a divergence- free vector such that Vq x s^^^^ = 0. 

Then, to see if the effect of transients can be suppressed by considering higher 
order terms, we continue to analyze the solution valid up to third-order Lagrangian 
perturbation theory (3LPT). Inserting the expansion Eq. ([9]) into Lagrangian Poisson 
equation, and using the results of the linear order and the second-order, we can obtain 
the third-order expression. It is more convenient to split the third-order displacement 
S^^^ into two parts as follows: 

S^'\q) = Sf\q) + Sf\q) , (27) 

where S^^ is from the cubic interaction among the linear perturbations and S^"* is from 
the interaction between the linear and the second-order perturbations. 
Then, the part for S^'^ is constrained by 



Fa{T) - a{T)Fa{r) = -2a{T)D{rf , (28 
^(1) 



V, ■ = det(.« ) . (29) 
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In the E-dS model, by substituting a(r) = 6r^^ into Eq. (1251) . the growing mode 
solution can be obtained as 

In addition to this, the solutions satisfying Eq. without the source term like 

Fa-(t) =C3t'/' + C4rS (31) 

where C3 and C4 are constants, also satisfy Eq. fl28l) and they serve as transients, too, 

unless the initial condition are set appropriately. 

However, compared with the growing mode, these transients are suppressed by 
and they are less problematic than the transients related to the second-order 

perturbation, which are suppressed by a~^. 

From Eq. fl29l) we can obtain the spatial part of Sa as 

where Ra is a divergence- free vector such that Vq x Sa = 0. 

(s) 

On the other hand, the part is constrained by 

F,(r) - a(r)F,(r) = -2«(r)Z}(r)[^(r) - ^(r)^] , (33) 

.(1) J2)" 



(V,-.«)(V,-.(^))-.S.gi . (34) 



(3) 

From the same discussion for the part Sa , we can acquire the growing mode solution 
in the E-dS model, 

as well as the spatial part of Sjj^\ 

) = \ (V, ■ .(^)) - (s(^) . V,) .(^) 

+ Rf , (36) 

f3) f3^ 

where Rl is again a divergence- free vector such that Vg x = 0. 
3. Statistics 

Here, we introduce the statistical quantities of the density fields on which transients 
provide the impact at large scales. For this purpose, a one-point probability distribution 
function of the density fluctuation field P{S) (PDF of the density fluctuation) which 
denotes the probability of obtaining the value 6 plays an important role. If 5 is a 
random Gaussian field, the PDF of the density fluctuation is determined as 
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where cr^ = <^(5 — (5))^) is the dispersion and ( ) denotes the spatial average. 

Since hnear growing modes are reproduced almost exactly by the ZA initial 
condition, the power spectrum, which is the quantity obtained by linear perturbation 
theory, is not affected by transients. Therefore, the expected tools to detect transients 
at large scales are the cumulants of the one-point PDF of the density fluctuation field 
whose orders are higher than two which become nonzero for the distribution deviating 
from Gaussian. In the structure formation, non-Gaussianity are generated because of the 
non-linear dynamics of the fluctuations, even though 5 is initially treated as a random 
Gaussian field, as a result of the generic prediction of inflationary scenario. 

For this purpose, we concentrate on the third and fourth order cumulants, which 
are defined as < 5^ >c=< >, < >c=< > — 3cr^ display asymmetry and non- 
Gaussian degree of "peakiness", respectively, for a given dispersion [H [2]. Since it is 
known that the scaling < 5" >cOC o"^""^ holds for weakly non-linear regions during the 
gravitational clustering ^llj from Gaussian initial conditions, we introduce the following 
normalized higher-order statistical quantities : 



skewness : 7 



kurtosis : ri = — ^ . 

The merit of adopting these definitions is, as stated above, that they are constants in 
weakly nonlinear stage which are given by Eulerian linear and second-order perturbation 
theory P, [TT]. For example, in the E-dS model smoothed with a spherical top-hat 
window function, 

^^Si^nx-xcosxl^ (38) 
the skewness and the kurtosis are given by 

7 = y + Z/i + 0{a') , (39) 
60712 62 7 o 2 2, , , 



where 



1323 3 " 3"' 3' 



with smoothing scale R |11] . 

For this form of the skewness and the kurtosis, the effects of transients at large 
scales from the ZA initial condition is also investigated by [161 lEj as 

6 12 

7tran = ~ 77 + 7777779 ' (^2) 



5a 35a'^/^ ' 
816 28?/i 184 1312 



35a 5a 75a2 245aV2 

8?/i 1504 192 , , 

H — \ f43) 

5a7/2 4725a9/2 ^ 1225a7 ' ^ ' 
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For the initial condition, instead of the ZA, if we adopt the one based on 2LPT, 
transients in the skewness (1421) is expected to be vanished at large scales, since 2LPT 
can provide appropriate initial condition up to second order related with the skewness. 
Similarly, for the one based on 3LPT, transients in the kurtosis (1431) is also expected to 
be vanished at large scales. 

Essentially, in the following section, we examine the effects of transients in the 
ACDM model relying on A^-body simulation from the initial conditions based on 2LPT 
and 3LPT. Since we investigate numerically, we can go into the highly nonlinear region, 
in which the analytic estimates of transients like Eqs. (142|) and (143!) can not be performed. 



4. Numerical Calculations 



In this section, we calculate the statistical quantities introduced in the previous section 
in the ACDM model based on A^-body simulations. For setting up the initial conditions, 
we use COSMICS code [29] which generates primordial Gaussian density field usually 
based on the ZA. We consider the case with those based on 2LPT and 3LPT, too. 
COSMICS package consists of 4 applications. GRAFIC generated Gaussian random 
density fields (density, velocity, and particle displacements) on a lattice. Both the 
velocity and the displacements are related to each other. 

GRAFIC automatically selects the output redshift by the maximum density 
fluctuation on a grid 5max for a given set of cosmological parameters. In order to obtain 
the initial redshift, we adopt the following cosmological parameters at the present time 
[z = 0) which are given by WMAP 3- year result [30] : 

1]„ = 0.28, (44) 

= 0.72, (45) 

Ho = 73 [km/s/Mpc] , (46) 

as = 0.74. (47) 

The averaged relation between the input maximum density fluctuation and the output 
redshift is shown in Table. [H The initial redshift is set by the input maximum density 
fluctuation. Because we set random Gaussian fluctuation, the initial redshift is not 
fixed. 

From the initial conditions set up above, we follow the evolution of the particles 
based on N-body simulation. The numerical algorithm is applied by particle-particle 
particle-mesh (P^M) method [6] which was developed by Gelb and Bertschinger. The 
numerical code we use is written by Bertschinger. For N-body simulations, we set the 
parameters as follows: 

Number of particles : A^ = 128^ , 

Box size : L = 128/i"^Mpc (at ^ = 0) , 

Softening length : e = 50/i^^kpc (at z = 0) . 



Transients from initial conditions based on Lagrangian perturbation theoryin N-body simulationslO 



Table 1. The averaged redshift and its dispersion of the initial condition selected by 
the GRAFIC code for a given maximum density fluctuation. 



^max 




z 




1.0 


22, 


.821 


0.868 


0.5 


32 


.833 


1.670 


0.2 


83, 


.806 


4.256 



For the simulations, we use 50 samples for an initial condition. After the 
calculations, in order to avoid the divergence of the density fluctuation in the limit 
of large k, however, just for a technical reason, it is necessary to consider the density 
field pm{x; R) at the position x smoothed over the scale R, which is related to the 
unsmoothed density field Pm{x) as 

pU^; R)= J d^yW{\x - y\-R)p^{y) , (48) 

where we use the top-hat spherical window function by Eq. (l38l) . Throughout this paper, 
we choose the smoothing scale R = 2/i~^Mpc (at z = 0). 

Crocce, Pueblas, and Scoccimarro [18j analyzed the non-Gaussianity of the density 
field with both the skewness and the kurtosis. They changed the smoothing scale R and 
compared the non-Gaussianity with two different initial conditions which are based on 
the ZA and 2LPT, respectively. Then they showed that the difference of the skewness 
and the kurtosis for R = 2h~^ Mpc becomes several percent at 2; = 0, when the initial 
condition is imposed at z = 24. 

Here, in addition to the ZA and 2LPT, we also analyze the non-Gaussianity with 
the initial conditions based on 3LPT, by slightly modifying the COSMICS so that 
initial conditions are set up by not only the ZA but also 2LPT and 3LPT. Using 
output files which describe the ZA displacement from COSMICS, we compute 2LPT 
and 3LPT displacements and velocities using Eqs. ( l26ll . (1321) . and ( l36ll . Here we set 

= Rf = i2g) = 0. The vectors R^'\Rf, and R^ were introduced for rotation- 
free condition. Because we set these vectors as zero, the rotation-free condition slightly 
violates. 

The main purpose is to see the impacts of transients from 2LPT initial conditions 
since 3LPT initial conditions are expected to provide exact results for the kurtosis. We 
investigate this also by changing the initial time dependence which corresponds to the 
initial maximum density fluctuation for a given Fourier mode. We expect this can clarify 
the range where N-body simulations from 2LPT initial conditions are accurate enough. 

First, we analyze the case in which the initial condition is imposed after the density 
fluctuation becomes relatively large (Smaxiti) = 1). 

Before analyzing the non-Gaussianity, we consider the evolution of the dispersion 
of the density distribution (Fig. [1]). We can see that the ZA is a good approximation 
initially, but it gradually deviates from the one from 2LPT initial condition because 
of the nonlinear effects. The final difference between the dispersions with the initial 
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Figure 1. The dispersion of tlie density distribution from N-body simulation {R = 
2h~^ Mpc) with different initial conditions. The difference between dispersion with 
the initial conditions based on the ZA and 3LPT is less than 4% at z = 0. 



conditions based on the ZA and 2LPT is less than 4%. There is almost no difference 
between the dispersions in density with the initial conditions based on 2LPT and 3LPT. 
Obviously, this confirms the strength of ZA, that can reproduce the growth of the density 
fluctuations well within the quasi-nonlinear regime. 

However we also see that when it comes to the non-Gaussianity which include the 
skewness and higher-order moments, the validity of the ZA is very limited (Fig. [2]). The 
difference between the skewness and the kurtosis with the initial conditions based on 
the ZA and 3LPT become about as much as 15% and 40% (at z = 2), respectively. 
Furthermore, the difference between them based on 2LPT and 3LPT becomes also 
significant. This means that the simulations with the initial condition based on 2LPT 
does not provide the correct skewness because nonlinearity at this initial time is 
significant. Without doubt, in such a case in which the initial condition is imposed 
at ^ ~ 23 like [1^, it is not enough to consider the initial condition based on 2LPT and 
we can acquire much more accurate values with the initial condition based on 3LPT. 

Next, we analyze the cases in which initial conditions are imposed for smaller density 
fluctuations. We reduce the maximum value of the density fluctuation Smaxiti) at initial 
time to 0.5, and 0.2. Decreasing 6max will lead to less analytically evolved solutions 
and more numerical evolution and hence is expected to decrease the effect of transients 
which are due to analytical evolution. For the case Smaxiti) = 0.5, we can see that 
the difference between the non-Gaussianity with the initial conditions based on the ZA 
and 3LPT is much smaller than the case Smaxiti) = 1.0 (Fig. [3]). This is because for 
the smaller density fluctuation, the higher order effect in LPT becomes less efficient 
and the ZA is better approximation for the initial condition. We can also see that the 
non-Gaussianity with the initial condition based on the ZA approaches those with the 
initial condition based on 2LPT and 3LPT at late time. This means that transients 
disappear until that. 

It is worth noting here is that for Smaxiti) = 0.5, there is no difference between the 
skewness with the initial conditions based on 2LPT and 3LPT since even 2LPT provides 
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Figure 2. The non-Gaussianity of the density distribution from N-body simulation 
{6„iax{ti) = 1.0, i? ~ 2h^^ Mpc) with different initial conditions, (a) The skewness 
of the density distribution. The difference between the skewness with the initial 
conditions based on the ZA and 3LPT is less than 14%. Furthermore, the skewness 
with the initial condition based on 2LPT does not coincide with that based on 3LPT 
obviously because of the high nonlinearity of the initial time, (b) The kurtosis of the 
density distribution. The difference between the kurtosis with the initial conditions 
based on ZA and 3LPT is less than 40%. Furthermore, the kurtosis with the initial 
condition based on 2LPT does not coincide with that based on 3LPT obviously. 

accurate skewness for this initial time with small nonlinearity. For the kurtosis, the 
difference with the initial conditions based on 2LPT and 3LPT is also very small ( at 
most less than 2%), because transients with 2LPT decreases faster than that with the 
ZA. This difference disappears until late time as the transients vanishes at that time. 

Furthermore, we examine the case Smaxiti) = 0.2 (Figs. H]). Both of the cases, the 
difference between the non-Gaussianity with the initial conditions based on 2LPT and 
3LPT almost disappears, as is expected from the results of the case Smaxiti) = 0.5. 
Even for the kurtosis, we cannot see the difference because before z = 5, transients with 
2LPT completely disappear. 

What is quantitatively important here is that the difference between the skewness 
and the kurtosis of the density distribution with the initial conditions based on the ZA 
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Figure 3. The non-Gaussianity of the density distribution from N-body simulation 
{5,riax{ti) — 0.5, i? — 2h^^ Mpc) with different initial conditions, we also show the 
error bars, (a) The skewness of the density distribution. The difference between the 
skewness with the initial conditions based on the ZA and 3LPT is less than 10%. 
Because of the suppressions of transients at late time, this difference disappears at 
late time. Notice that there is no difference between the skewness with the initial 
conditions based on 2LPT and 3LPT since 2LPT provides almost exact skewness for 
this initial time with small nonlinearity. (b) The kurtosis of the density distribution. 
The difference between the kurtosis with the initial conditions based on the ZA and 
3LPT is less than 25%. Because of the suppressions of transients at late time, this 
difference disappears at late time. Notice that the difference between the kurtosis with 
the initial conditions based on 2LPT and 3LPT is less than 2% and almost disappears 
at late time because transients with 2LPT decrease much faster than the one with the 
ZA. 



and 2LPT still survives for such a small value of Smax{ti) = 0.2 until z ~ 1. This means 
that it takes a long time for transients with the ZA to completely disappear. 

By now, we only show the results after z = 5, even though we have actually done 
the numerical calculations from about z = 20. For the reader who are interested in the 
thorough evolution of the statistical quantities, we show one example with 6maxiti) = 0.5 
(Fig. 5). 

To finish this section, it is worth commenting on the error bars in the curves (Fig[3]). 
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(b) 



Figure 4. The non-Gaussianity of the density distribution from N-body simulation 
{^max{ti) = 0.2, i? — 2h^^ Mpc) with different initial conditions, (a) The skewness of 
the density distribution. There is no difference between the skewness with the initial 
conditions based on 2LPT and 3LPT since 2LPT provides almost exact skewness for 
this initial time with small nonlincarity. The difference between the skewness with the 
initial conditions based on the ZA and 2LPT is less than 4% but remains until z 1 
because of transients, (b) The kurtosis of the density distribution. The difference 
between the kurtosis with the initial conditions based on 2LPT and 3LPT become 
negligible because transients with 2LPT have disappeared until 2^5. The difference 
between the kurtosiss with the initial conditions based on the ZA and 2LPT is less 
than 10% but remains until z ^ 1 because of transients. 

Although we only show error bars the Smaxiti) = 0.5 case, their values are almost 
independent of the initial conditions (ZA, 2LPT or 3LPT), hence hold for the other 
simulations. In Table [2] and Table [3l we show the concrete values for the errors. 

5. Nonlinear effects in top-hat spherical symmetric model 

According to the numerical calculation in the previous section, we see that even though 
the difference between the non-Gaussianity with the initial conditions based on the ZA 
and 2LPT can not be removed well by taking the initial time at sufficiently high-^ era, 
it can be significantly decreased between those based on 2LPT and 3LPT. Especially, 
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Figure 5. The non-Gaussianity of the density distribution from N-body simulation 
{5raax{ti) — 0.5, i? — 2h^^ Mpc) with different initial conditions. They are based on 
the same calculation as in Fig. [Sj but show the evolution of the statistical quantities 
from z = 20. (a) The skewness of the density distribution, (b) The kurtosis of the 
density distribution. 

Table 2. The error of the statistical quantities at z = 2 from N-body simulation 
{Smax{ti) = 0.5, i? = 2h^^ Mpc) with different initial conditions. 



Initial condition 


Dispersion 


Skewness 


Kurtosis 


ZA 


3.2% 


19.6% 


74.1% 


2LPT 


3.4% 


21.8% 


82.6% 


3LPT 


3.5% 


22.1% 


84.3% 



the difference between the kurtosis which initially includes transients with 2LPT has 
disappeared until 2; ~ 2, if we set the initial time for A^-body simulation around z ~ 33 
(Table p. 

For the interpretation and further analytical justification for the results, in this 
section, we reconsider a simpler situation modeled by a top-hat spherical symmetric 
collapse with constant density for which both the exact solution and Lagrangian 
perturbative expansion are obtained analytically [121 [25]. Another motivation 
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Table 3. The error of the statistical quantities at z = from N-body simulation 
{Smax{ti) = 0.5, R = 2h~^ Mpc) with different initial conditions. 



Initial condition 


Dispersion 


Skewness 


Kurtosis 


ZA 


6.9% 


17.8% 


46.6% 


2LPT 


7.0% 


17.5% 


45.1% 


3LPT 


7.0% 


17.5% 


45.2% 



for the symmetric collapse, is that it provides most severe constraints for Lagrangian 
perturbation theory, in general, as shown by Yoshisato et al. [31] . 

In this section, we reset the definition of the scale factor a. If such a shell is in the 
Einstein-de Sitter universe, the equation of motion for a spherical shell is given by 



d 



— a 



dt 



^dx\ 



2a X 



X / 



(49) 



where a{t) oc t"^^^ is a scale factor, a; is a comoving radial coordinate and Xq = x(to) 
For the derivation of Eq. P9|) . see around Eq. (17) in [2T] . 

Under the initial condition 5 = a for a — 0, Eq. (H9l) can be integrated as. 



\ da 



1 

R 



(50) 



where R{t) = a{t)x/xQ is a physical particle trajectory. It is known that the exact 
solution for the spherical collapse (Eq. (l50!l ) can be parametrized as follows: 



me) = ^(1 



a{9) 



cos I 



— sm I 



2/3 



(51) 
(52) 



From Eqs. (15T!) and (!52|) . we can also obtain the density contrast given hj S = {xq/x)^ 
as [El EH [25], 

9(6 -sine? 



6{x) 



1. 



(53) 



2(1 - cos^) 

Together with Eqs. (!52l) . (l53l) and the definition of the density contrast, we can express 
R{t) exactly in terms of the function a{t). 

Obviously, this exact expression can be expanded in terms of a like 



Rit) = Ro 



1 - J2i-^)'Cka' 



k=l 



(54) 



where Ck are perturbative coefficients. It can be shown that the first three terms of this 
expansion are given as 



23 
1701 



(55) 



Since Eq. (IMIl is an expansion with respect to the displacement from the 
homogeneous distribution, this is nothing but the Lagrangian perturbation. 
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On the other hand, for a general symmetric spherical symmetric collapse, the 
density fluctuation in Eulerian picture at r can be related with the perturbation variable 
in Lagrangian picture as follows: 

3{r) =(l- ^s{r)\ - 1, (56) 



where s(r) means the radial component of Lagrangian displacement which can be 
expanded as 

s(r) = s« W + + s(=^)(r) + 0{e^) . (57) 

Compared with Eq. (IMl) and Eq. together with the definition of the density 
contrast, we can obtain the following relations: 

k=l 

and applying the expansion given by Eq. fl57j) this yields 

|-.«(r) = ia, (59) 

= La' = U |-s«(r)V , (60) 



dr 21 7 \dr 



9r ' ' 1701 63 

This spherical model forms caustics at a = 3 in the ZA. Here we note that this is an ideal 
model, and the value of a which forms caustics is different from the actual Universe. 

In GRAFIC code, the initial conditions are automatically set when the maximum 
density fluctuation at any lattice point becomes greater than some value. If the point is 
approximated as the center of the spherical collapse, we can derive the relation between 
the density fluctuation and the linear Lagrangian displacements from Eq. fl56l) . 

^s^^\r) = l-{l + 6{r)y^/\ (62) 
or 

The linear Lagrangian displacements at this time becomes 
d 

^ 0.206299 for 6 = 1, (63) 

or 
d 

—s^^\r) ~ 0.126420 for 6 = 0.5, (64) 
Or 

d 

—s^^\r) ~ 0.058964 for 6 = 0.2. (65) 
or 

Therefore based on the discussions in this section, we can calculate the effect of 2LPT 
and 3LPT at initial time (Table H]). 

Here we regard the actual Universe with the above result. Clearly, by regarding the 
top-hat spherical collapse considered in this section as the generation of the overdense 
regions through gravitational instability in the Universe, this analysis provides another 
explanation for the results in the previous section. For example, for the initial condition 
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Table 4. The effect of tlie second- and tliird-order perturbation for the density 
fluctuation evaluated at a given 5 during the top-hat spherically symmetric collapse. 
5*-*' denotes that this quantity is calculated by i—th order Lagrangian perturbation 



theory. 


6 




5(3)/(5(l)+5(2)+5(3)) 


1.0 


14.4% 


2.34% 


0.5 


7.17% 


0.73% 


0.2 


2.86% 


0.14% 



^max = 1, the effect of 3LPT lias become significant, while for 5max = 0.5 and 6max = 0.2 
it is only below 1% and eventually negligible. On the other hand, the effect of 2LPT is 
still significant even we choose the initial condition with smaller density fluctuation like 

^max 0.2. 

6. Summary 

Recently, observations have improved the precision of cosmological constraints 
tremendously. In order to compare theoretical predictions with observational results, 
we should carry out detailed numerical simulations with high accuracy. In this context, 
in order to follow the nonlinear evolution, one needs to consider carefully how to set 
initial conditions of A^-body simulations very carefully. In the standard method, the 
initial conditions are set up by the Zel'dovich approximation (ZA) in the quasi-nonlinear 
regime. The ZA has the nice property that by using linear perturbation theory allows 
one to describe accurately the quasi-nonlinear density field, in the case of small non- 
linearity. 

On the other hand, when constraining cosmological models, non-Gaussianity of the 
density fluctuations generated by the nonlinear dynamics will become important. It is 
well known that the ZA is not sufficient to measure non-Gaussianity from higher-order 
statistics, for example, the skewness and kurtosis. The reason is the existence of so- 
called transients, that is, the most dominant decaying mode arising from our ignorance 
of the initial conditions. Until the effects of transients have disappeared, we can not 
accurately reproduce higher-order quantities like the skewness and the kurtosis. 

Actually, recently, Crocce, Pueblas, and Scoccimarro analyzed the impact 
of transient on the second- and higher-order statistics quantitatively for A^-body 
simulations by comparing the initial conditions based on the ZA and second-order 
Lagrangian perturbation theory (2LPT). They show that even though there are also 
transients from 2LPT initial conditions, the effects of transients are suppressed compared 
to the ZA case. 

In this paper, as a natural extension of [18], in addition to the ZA and 2LPT, 
we also analyze the non-Gaussianity with the initial conditions based on third-order 
Lagrangian perturbation theory (3LPT). Since 3LPT initial conditions are expected to 
provide exact results for the kurtosis in the weakly nonlinear region, we can evaluate 
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the impact of transients from 2LPT initial conditions. The goal is to clarify the range 
where A^-body simulations from 2LPT initial conditions are accurate. 

When we set the initial condition of the maximal density contrast reaching unity 
i^max ~ 1) at 2; ~ 20, differences in the predicted non-Gaussianity between the three 
different initial conditions (ZA, 2LPT and 3LPT) are readily apparent. Since the 
accurate value of the skewness is reproduced by 2LPT for the weakly nonlinear regime, 
this disagreement is because of the nonlinearity at this initial time. Therefore, 2LPT 
initial conditions are not sufficient to set up the initial condition for A^-body simulations. 
There is also no guarantee that 3LPT initial conditions are accurate enough for precise 
determination of cosmological parameters. 

Next, to avoid the problem with the initial nonlinearity, we set the initial conditions 
at z ~ 30, corresponding to the maximal density contrast is about half, {6max ~ 0.5). 
In this case, as is expected from the predictions of 2LPT in the weakly nonlinear region, 
we find no difference between the 2LPT and 3LPT results, confirming that the 2LPT 
initial conditions produce accurate skewness. For the kurtosis, the difference between 
the initial conditions with 2LPT and 3LPT is also very small (2%). On the other hand, 
the difference of non-Gaussianity with initial conditions based on the ZA and 3LPT is 
large (10% for the skewness and 25% for the kurtosis) until 2; ~ 1. This shows that 
transients from initial conditions with 2LPT have less impact than the ones with the 
ZA initial conditions. 

This tendency is also obtained for the initial conditions set at z ~ 80 when the 
maximal density contrast is about two-tenths, {6max ~ 0.2). From the numerical 
calculations for the kurtosis, we can see that the effects of transients from 2LPT initial 
condition have completely disappeared until z ~ 5. However, there is still significant 
differences in the predicted non-Gaussianity with initial conditions based on the ZA and 
3LPT (4% for the skewness and 10% for the kurtosis) until z ~ 1. 

Therefore, for sufficiently early z > 30 initial conditions, the effects of transients 
on kurtosis from 2LPT initial conditions become negligible until roughly 2; ~ 3, while 
those from the ZA initial conditions survive until z ~ 1. As long as one considers typical 
A^-body simulations, which start at 2; ~ 50, the predicted statistics are accurate enough 
up to the forth-order (kurtosis) using 2LPT initial conditions. 

We further investigated our above results semi-analytically by considering the 
simpler situation modeled by a top-hat spherically symmetric collapse with constant 
density. We can show that in this situation, the difference between the impact of 2LPT 
and 3LPT initial conditions becomes almost negligible for the initial maximum density 
contrast much less than the half {6max ~ 0.5), while the difference between the impact 
of the ZA and 2LPT initial conditions is still greater than 2% for the initial maximum 
density contrast about two-tenth {6max ~ 0.2). 

Finally, we briefly mention the observational consequences of our results. It is 
known that weak lensing surveys can potentially provide us with precision maps of the 
projected density up to redshifts around 1 [32l [33l [Ml [351 EE]. Even though we need 
another step in obtaining the convergence field which can be written as the projection of 
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the matter density along the hne of sight, the skewness and kurtosis of the convergence 
field can be tested via weak lensing surveys. Our comparison of the initial conditions 
should prove useful extracting cosmological parameters from observational data. 
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